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ABSTRACT 

An  analytical  approach  to  the  problem  of  replenishing  task  forces 
at  sea  is  investigated,  using  a  random  walk  model.   The  type  of 
replenishment  operations  considered  consists  of  two  replenishment 
ships  of  different  types,  each  replenishing  L  combatant  ships  of 
different  types.   The  combatant  ships  are  initially  distributed  into 
queues  of  length  M  and  N  between  the  two  replenishment  ships.   As  each 
combatant  ship  completes  its  first  replenishment,  it  then  enters  the 
queue  of  the  other  replenishment  ship.   Thus  each  combatant  ship  is 
replenished  by  both  replenishment  ships.   Replenishment  times  for  the 
combatant  ships  are  assumed  to  be  independent,  positive  valued  random 
variables  whose  distributions  are  dependent  upon  both  replenishment 
ship  type  and  combatant  ship  type.   The  replenishment  of  the  task 
force  is  completed  whenever  2L  replenishments,  two  for  each  combatant 
ship,  are  completed.   The  total  replenishment  time  is  shown  to  be  the 
maximum  of  certain  partial  sums  of  the  individual  replenishment  times. 
A  general  expression  for  the  distribution  function  of  the  total  re- 
plenishment time  is  sought  by  various  methods  of  analysis.   A  computer 
simulation  is  used  to  gain  a  better  understanding  of  the  form  of  the 
distribution  function  of  total  replenishment  time  for  a  specific  under- 
way replenishment  operation.   Some  examples  are  discussed. 
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I.   THE  UNDERWAY  REPLENISHMENT  OPERATION 

The  underway  replenishment  operation  is  the  means  by  which  needed 
supplies  are  provided  to  Navy  fleet  units  while  underway  at  sea.   It 
provides  these  units  with  the  capability  of  operating  at  great  distances 
from  Navy  shore  facilities,  enables  them  to  remain  on  station  for  ex- 
tended periods,  extends  their  period  of  effective  combat  readiness, 
and  limits  diversion  from  their  assigned  mission.   The  supplies  pro- 
vided may  include,  but  need  not  be  limited  to,  petroleum  products, 
missiles,  ammunitions,  stores,  and  provisions.   However,  no  matter 
what  is  to  be  transferred  between  the  transferring  replenishment  ships 
and  the  receiving  combatant  ships,  several  points  are  apparent,  the 
underway  replenishment  operation  increases  the  vulnerability  of  the 
engaged  units,  restricts  their  mobility,  and  diverts  the  involved 
combatant  ships  from  their  principal  missions.   The  Navy  is  continually 
striving  for  the  development  of  techniques  and  methods  which  will 
reduce  the  limitations  imposed  on  fleet  operations  by  the  necessity 
for  underway  replenishment.   Principally,  the  Navy  desires  to  reduce 
the  vulnerability  of  the  engaged  units,  increase  the  unit  mobility 
while  engaged,  reduce  the  time  required  to  complete  replenishment, 
and  limit  interference  with  the  assigned  missions  of  the  combatant 
units . 

To  attain  these  goals  many  improvements  of  a  material  nature  have 
been  accomplished.   The  introduction  of  faster  replenishment  ships, 
such  as  the  combat  stores  ship  (AFS),  reduces  the  vulnerability  and 
increases  the  mobility  of  the  engaged  units  by  allowing  underway 
replenishments  to  be  conducted  at  higher  speeds  than  before.   The 


introduction  of  mu it i -commodity  replenishment  ships,  such  as  the  fast 
combat  support  ship  (AOE)  and  the  fleet  replenishment  tanker  (AOR) , 
reduces  the  total  number  of  replenishments  required  since  several  needs 
can  be  filled  during  a  single  replenishment.   The  introduction  of  the 
helicopter  as  a  transfer  agent  affords  greater  separation  of  task  force 
units,  reduces  the  demands  for  alongside  replenishments,  and  minimizes 
interference  with  the  assigned  mission  of  the  involved  combatant  ship. 
The  introduction  of  new  wire-highline  transfer  systems,  such  as  the 
fast  automatic  schuttle  transfer  cargo  handling  system  (FAST),  provides 
a  higher  rate  of  cargo  transfer  and  an  increased  load  capacity,  thereby 
reducing  the  required  time  alongside.   In  addition,  internal  cargo 
handling  improvements  in  both  the  receiving  and  transferring  ships 
have  resulted  in  reduced  breakout  and  strikedown  times,  thus  again 
resulting  in  reduced  replenishment  times.   All  these  improvements  have 
had  one  element  in  common.   Their  accomplishment  has  resulted  in  a 
reduction  in  the  actual  replenishment  time  for  an  individual  ship,  and 
a  subsequent  reduction  in  the  total  replenishment  time  for  the  involved 
task  force . 

Improvements  have  not  been  limited  to  those  of  a  material  nature. 
The  underway  replenishment  operation  has  also  been  viewed  as  a  problem 
in  operations  analysis.   The  objective  of  such  an  approach  is  a  better 
understanding  of  the  underway  replenishment  operation  as  viewed  in  an 
analytical  manner. 

The  underway  replenishment  operation  can  be  viewed  in  many  analytical 
contexts,  but  more  readily  lends  itself  to  representation  as  a  multi- 
server  queueing  model  with  a  finite  number  of  customers.   The  number  of 
servers  can  be  established  by  the  number  of  replenishment  ships  involved, 
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while  the  number  of  customers  can  be  established  by  the  number  of  ships 
in  the  task  force  scheduled  for  underway  replenishment  operations.   It 
should  be  noted  that,  under  normal  conditions,  a  definitive  sequence  of 
individual  replenishments  has  been  established  prior  to  the  commencement 
of  the  underway  replenishment  operation.   This  is  done  as  a  method  of 
optimizing  replenishment  ship  utilization,  while  hopefully  minimizing 
total  replenishment  time  for  the  task  force.   Thus,  while  the  underway 
replenishment  operation  may  be  viewed  in  terms  of  a  queueing  problem, 
it  has  definite  job-shop  scheduling  problem  overtones. 

One  aspect  of  a  typical  queueing  model  must  be  ignored  for  the 
underway  replenishment  operation.   It  is  normal  to  consider  all  service 
times  to  have  come  from  the  same  probability  distribution  in  a  typical 
queueing  model.   Since  many  different  types  of  replenishment  and  com- 
batant ships  are  involved  in  a  typical  underway  replenishment  operation, 
this  seems  to  be  an  unreasonable  assumption  for  any  realistic  model  of 
the  operation.   Analysis  of  actual  underway  replenishment  data  has 
supported  this  contention  and  has  yielded  approximations  for  the 
appropriate  distribution  functions  for  various  typical  underway 
replenishment  operations. 
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II.   THE  UNDERWAY  REPLENISHMENT  MODEL 

The  basic  underway  replenishment  model  that  was  developed  by 
Waggoner  [1]  and  later  expanded  by  Patterson  [2]  is  the  domain  of 
inquiry  for  this  thesis.   This  model  assumes  that  two  replenishment 
ships  and  L  combatant  ships  are  scheduled  to  conduct  underway  replenish- 
ment operations;  a  total  of  2L  replenishments  are  to  be  conducted,  L  by 
each  replenishment  ship.   Prior  to  the  commencement  of  the  replenishment 
operation,  the  L  combatant  ships  are  given  a  specific  replenishment 
sequence  and  assigned  to  queues  of  length  M  and  N,  where  M  +  N  =  L, 
associated  with  the  first  and  second  replenishment  ships,  respectively. 
For  each  combatant  ship  this  means  that  a  definitive  order  of  replenish- 
ment has  been  established  and  that  each  combatant  ship  has  initially 
been  assigned  to  one  of  the  two  queues.   In  addition,  this  means  that 
the  individual  combatant  ship  both  succeeds  and  precedes  the  same  ships 
in  each  queue  and  that  this  queueing  sequence  is  inalterable. 

The  movement  of  the  combatant  ships  during  the  underway  replenishment 
operation  is  envisioned  as  a  sequence  of  steps  in  a  random  walk.   This 
random  walk  is  performed  on  an  L  x  L  grid  on  which  the  coordinates  of 
each  grid  node  represent  a  possible  combination  of  the  numbers  of  com- 
batant ships  replenished  by  each  replenishment  ship,  i.e.,  the  node 
(i,j)  on  the  grid  indicates  that  i  combatants  have  been  replenished  by 
the  first  replenishment  ship  and  j  by  the  second  replenishment  ship. 
It  is  obvious  from  the  replenishment  sequence  that  certain  grid  nodes 
represent  infeasible  combinations;  it  is  infeasible,  for  example,  for 
the  first  replenishment  ship  to  have  completed  M  +  k  replenishments 
unless  the  second  replenishment  ship  has  completed  at  least  k,  where 
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k  =  1,2,..., N.   This  infeasibility  can  be  stated  in  constraint  equations 
If  the  abscissa,  Z, ,  of  the  grid  denotes  the  number  of  replenishments 
completed  by  the  first  replenishment  ship,  while  the  ordinate,  Z~ , of  the 
grid  denotes  the  number  of  replenishments  completed  by  the  second  re- 
plenishment ship.   The  following  equations  constitute  the  necessary 
constraints : 

Z  =  0 

Zl   =  M  +  Z2 

Zl   =  L 

z2  =  o   • 

Z2  =  N  +  Zx 
Z2  =  L 

All  feasible  nodes  either  lie  on  or  within  the  boundaries  established 
by  the  intersection  of  these  constraint  equations.   Figure  1  depicts 
an  L  x  L  grid  with  the  constraint  equations  superimposed. 

In  order  to  establish  a  common  terminology  for  later  discussion 
and  investigation  of  the  model,  various  definitions  and  theorems 
stated  by  Patterson  [2],  are  presented. 
Definition  1 

The  CONSTRAINT  BOUNDARY  of  the  random  walk  grid  is  defined  as  follows 
(see  Figure  1)  : 

a.  the  UPPER  BOUNDARY  consists  of  those  nodes  along  the  line 
Z~  =  L  with  abscissa  k,  where  M  <  k  <  L-l; 

b.  the  RIGHT  BOUNDARY  consists  of  those  nodes  along  the  line 
Z,  =  L  with  ordinate  h,  where  N  <  h  <  L-l; 

c.  the  LOWER  BOUNDARY  consists  of  those  nodes  along  the  line 
Z     =  M  +  Z2,  where  0  <  Z2  <  H-l;  and 
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Underway  Replenishment  Model 
Random  Walk  Grid 

FIGURE  1 
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d.   the  LEFT  BOUNDARY  consists  of  those  nodes  along  the  line 
Z2  =  N  +  Z  ,  where  0  <  Z  <  M-l. 
Definition  2 

A  STEP  is  a  horizontal  or  vertical  line  between  adjacent  nodes  of  the 
grid;  a  PATH  is  a  sequence  of  non-negative  steps  which  starts  at  the 
origin  and  terminates  at  (L,L)  and  does  not  cross  the  constraint 
boundary;  and  X.  and  Y.  denote  the  times  required  to  complete  the  ith 
horizontal  and  jth  vertical  step,  respectively,  of  any  path. 

From  Definition  2  it  is  evident  that  X.  is  the  ith  replenishment 
time  for  the  first  replenishment  ship,  while  Y.  is  the  jth  replenish- 
ment time  for  the  second  replenishment  ship. 
Definition  3 

The  set  of  all  feasible  paths  from  the  origin  to  (L,L)  is  decomposed 
into  four  mutually  exclusive  and  exhaustive  subsets: 

A:,  the  set  of  all  paths  that  pass  through  the  node  (L-1,L)  and 
do  not  pass  through  any  nodes  on  the  lower  boundary; 

B:   the  set  of  all  paths  that  pass  through  the  node  (L-1,L)  and 
pass  through  at  least  one  node  on  the  lower  boundary; 

C:   the  set  of  all  paths  that  pass  through  the  node  (L,L-1)  and 
do  not  pass  through  any  nodes  on  the  left  boundary;  and 

D:   the  set  of  all  paths  that  pass  through  the  node  (L,L-1)  and 
pass  through  at  least  one  node  on  the  left  boundary. 

In  addition,  subsets  B  and  D  can  be  further  decomposed: 

B,  :   the  subset  of  all  paths  in  B  which  last  pass  through  the 
lower  boundary  at  node  (M  +  k,k)  for  k  =  0,1,..., N-l 

D,  :   the  subset  of  all  paths  in  D  which  last  pass  through  the 
left  boundary  at  node  (h ,  N+h)  for  h  =  0,1,...,  M-l. 
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Theorem  1 

The  total  replenishment  time  is  the  actual  time  required  to  complete  a 

path  through  the  grid  and  is  equal  to  one  of  the  following  values: 

a.  the  total  length  of  time  required  for  any  path  in  A  is 

L 
A   i=l  X 

b.  the  total  length  of  time  required  for  any  path  in  B,  is 

L     k+1 
T   =    S   +  Z  Y.  for  k  =  0,1,... ,N-1; 
i=ttfk+l   j  =  l   L 

c.  the  total  length  of  time  required  for  any  path  in  C  is 

L 
T  =  S  Y  ;  or 
C   J-l   l 

d.  the  total  length  of  time  required  for  any  path  in  D,  is 

h+1         L 
T -  S  X.  +    S   Y.  for  h  =  0,1,. . . ,M-1. 
i=l  X   j=N+h+l  X 

Waggoner  and  Patterson  both  assumed  that  the  X.  and  Y.  were  inde- 
pendent, random  variables  from  the  family  of  exponential  probability 
distributions.   Waggoner  tacitly  assumed  that  the  replenishment  rate 
was  dependent  only  upon  the  replenishment  ship  type  and  that  replenish- 
ment rates  \    and  g.  could  be  used  for  all  replenishments  conducted  by 
the  first  replenishment  ship  and  the  second  replenishment  ship,  res- 
pectively.  Patterson  generalized  this  assumption  to  the  extent  that 
the  replenishment  rate  was  allowed  to  be  dependent  upon  both  the 
replenishment  ship  and  the  combatant  ship  types  involved.   Thus  re- 
plenishment rates  \.    and  \i.   were  assumed  for  the  ith  replenishment  by 
the  first  replenishment  ship  and  the  jth  replenishment  by  the  second 
replenishment  ship,  respectively.   It  should  be  noted  here,  however, 
that  Patterson  developed  Theorem  1  without  referral  to  his  assumption 
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of  exponentiality.   Thus  Theorem  1  holds  without  any  assumption  being 
made  regarding  the  distribution  of  the  individual  replenishment  times 
An  interesting  comparison  can  be  made  regarding  the  investigative 
approaches  of  Waggoner  and  Patterson.   Waggoner  approached  the  model 
by  considering  the  first  step  taken  and  the  subsequent  sequence  of 
steps  of  shortest  duration  taken  during  the  random  walk.   Thus  the 
exponentiality  assumption,  by  means  of  the  memoriless  property, 
insured  that  the  duration  of  each  step  continued  to  be  exponentially 
distributed.   This  greatly  simplified  his  derivation  of  the  Laplace 
transform  for  the  distribution  function  of  the  total  replenishment 
time.   On  the  other  hand,  Patterson  reversed  the  sequence  and  con- 
sidered the  last  step  taken  and  the  subsequent  sequence  of  steps  of 
longest  duration  in  the  reverse  order.   In  this  case,  as  noted  above, 
the  exponentiality  assumption  was  not  applied.   Patterson,  then  using 
the  exponentiality  assumption,  was  able  to  develop  the  Laplace  trans- 
form for  the  distribution  function  of  total  replenishment  time  only 
for  the  case  of  two  combatant  ships  initially  distributed  one  to  each 
replenishment  ship. 
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III.   TOTAL  REPLENISHMENT  TIME 


The  previous  section  has  intimated  that  total  replenishment  time 

can  be  one  of  certain  feasible  values.   It  was,  in  fact,  stated  that 

total  replenishment  time  was  either  T  ,  or  T    (for  some  integer  k 

A       u 

where  0  <  k  <  N-l),  or  T  ,  or  T    (for  some  integer  h  where  0  <  h  <  M-l) 
In  this  section,  by  means  of  the  following  theorem  and  proof,  it  is 
shown  that  total  replenishment  time  is  actually  the  maximum  of  those 
feasible  values. 
Theorem  2 

The  total  replenishment  time,  i.e.,  the  time  required  to  move  from  node 
(0,0)  to  node  (L,L)  of  the  random  walk  grid,  is  the  maximum  of  all 
feasible  values.   Let  T  denote  total  replenishment  time,  then 

T  =  max  (TA;  Tfik  for  k  =  0,1,..., N-l;  Tc  ;  T^  for  h  =  0 , 1 ,  .  .  .  ,M-1)  . 

Proof 

a.   Assume  that  the  path  taken  was  from  A.   Thus  it  must  be  shown 

that  T  >  T  k  for  k  =  0,1,..., N-l;  T  >  T  ;  and  T  >  T  h  for 

h  =  0,1,.. . ,M-1. 

1)  Because  the  path  is  in  A,  it  must  not  touch  the  lower 

boundary.   Thus  it  must  reach  the  line  Z~  =  k  +  1  (see  Figure  1)  before 

it  reaches  the  line  Z.  =  M+  k  for  any  k  =  0,1,..., N-l.   This  implies 

that 

M+k      k+1 

E  X.  >  £  Y.  for  any  k  =  0,1,...,N-1. 
1=1   L   j=l   J 

L 
Now  by  adding  the  partial  sum    Z   X.  to  both  sides  of  the  above 

i=M+k+l  ^ 
k 
inequality,  it  is  seen  that  T  >  T   for  any  k  =  0,1,..., N-l. 

A     B 
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2)  Because  the  path  is  in  A,  it  must  reach  the  upper  boundary 

before  it  reaches  the  right  boundary.   This  implies  that 

L        L 

S  X.  >  E  Y.. 

i=l  L   j  =  l  J 

Therefore,  T  >  T_. 

3)  Let  the  path  in  A  first  touch  the  line  Z~  =  N  +  h  (see 
Figure  1)  at  the  node  (r,N+h)  where  h  <  r  <  L-l  for  any  h  =  0,1,..., M-l, 
Since  the  line  Z»  =  N  +  h  must  be  reached  vertically,  the  (r  +  l)th 
horizontal  step  must  have  begun  before  the  node  (r,N+h)  is  reached. 

The  (N-RH-l)th  vertical  step  is  either  just  beginning  (h+1  <  r  <  L-l) 
at  node  (r,Nfh)  or  will  begin  on  completion  of  the  (r  +  l)th  horizontal 
step  (  r  =  h).   Since  the  upper  boundary  must  be  reached  before  the 
right  boundary,  this  implies  that  either 


Z  X.  >    S   Y.  if  h+1  <  r  <  L-l,  or 
i=r+l  X   j=Nfh+l  J 

L  L 

£  X.  >    T,       Y.  if  r  -  h. 
i=r+2  X        j=N+h+l  J 

In  either  case,  the  above  inequalities  imply  that 


E  X.  >    E   Y.  for  any  h  =  0,1,..., M-l. 
i=h+2  L   j=N+h+l  J 
h+1 
Now  if  the  partial  sum  £  X.  is  added  to  both  sides  of  the  last 

1=1 1  h 

inequality,  it  can  be  seen  that  T  >  T   for  any  h  =  0,1,..., M-l. 

A  D 

b.   Assume  that  the  actual  path  taken  was  from  B,  where  k  is  fixed 

k         k     n 
and  0  <  k  <  N-l.   Thus  it  must  be  shown  that  T„   >  TA ;  T„   >  T_   for 
—   —  B     A   B     B 

n  jt  k,  0  <;  n  <:  N-l;  TBk  >  T  ;  and  TT,k  >  T  h  for  h  =  0,1,..., M-l. 

1)  Because  the  path  is  in  B,  ,  it  must  reach  the  line  Z.  =  M  +  k 
(see  Figure  1)  before  reaching  the  Z_  =  k  +  1.   Thus 
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k+1      M+k 

S  Y.  >  £  X.  . 
j-1  J    1-1  X 

L 

Now  by  adding  the  partial  sum        X.  to  both  sides  of  the  above 

i=M+k+l  X 

inequality,  it  can  be  seen  that  T„  >  T. . 

2)  Because  the  path  is  in  B,  ,  it  must  pass  through  the  node 

(M+k,k).   Thus  the  (M+k+l)th  horizontal  step  begins  at  the  same  time 

as  the  (k+2)th  vertical  step.   Since  the  path  cannot  pass  through  any 

node  (M+u,u)  where  k  <  u  <  N-l,  it  must  reach  the  line  Z_  =  u  +  1 

(see  Figure  1)  before  reaching  the  line  Z,  =  M+u .   These  conditions 

imply  that 

M+u        u+1 
E    i  >   E  Y.,  for  k  <  u  <  N-l. 

i=M+k+l     j=k+2  J 

L         k+1 

Now  by  adding  the  two  partial  sums    E   X.  and  £  Y.  to  both  sides 

i=M4-u+l  X  j=l   J 

k     u 
of  the  above  inequality,  it  can  be  seen  that  T   >  T   for  k  <  u  <  N-l. 

D  D 

3)  Let  the  path  in  B,  first  touch  the  line  Z„  =  v+1  (see 
Figure  1)  at  the  node  (s,v+l)  where  0  <  s  <  M+v  and  0  <  v  <  k.   The 

(v+2)th  vertical  step  begins  at  this  node,  while  the  (s+l)th  horizontal 
step  either  just  begins  (s  =  Mfv)  or  has  already  begun  (0  <  s  <  M+v-1). 
Since  it  is  necessary  that  the  line  Z„  =  k+1  (see  Figure  1)  not  be 
reached  until  after  the  line  Z,  =  M  +  k  is  reached,  this  implies  that 


either 


k+1  M+k 
£  Y.>   Z  X.  if  s  =  M  +  v  or 

j=v+2  J  i=s+l  X 

k+1  M+k 
£  Y.>   E  X.  if  0<  s  <M+v-l. 

j=v+2  J  i=s+2  L 


In  either  case,  this  implies  that 
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k+1  M+k 

£     Y.   >         £       X.    for   0  <  v  <  k. 
j=v+2   J        i=M+v+l   L 

v+1  L 

Now  if  the  partial  sums  S  Y.  and        X.  are  added  to  both  sides  of 

j=l  J     i=M+k+l  X 

k     v 
the  last  inequality,  it  can  be  seen  that  T   >  T   for  0  <  v  <  k. 

4)  Because  the  path  is  in  B,  ,  the  (k+2)th  vertical  step  and 
the  (M+k+l)th  horizontal  step  begin  at  the  same  time.   Since  the 
upper  boundary  must  be  reached  before  the  right  boundary,  this  implies 
that 

L  L 

£   X.  >   2  Y.. 
k=M+k+l  X        j=k+2  J 

k+1 

Now  by  adding  the  partial  sum  £  Y.  to  both  sides  of  the  above 

j  =  l   J 

inequality,  it  can  be  seen  that  T   >  T  . 

B      L 

5)  Because  the  path  is  in  B,  ,  it  cannot  reach  the  line  Z_  =  k+1 

(see  Figure  1)  before  reaching  the  line  Z,  =  M+k,  which  implies  that 

k+1      M+k 

£  Y.  >  I  X.  . 
j=l   J    i=l   1 

This  in  turn  implies  that 

k+1      h+1 
(1)   S  Y.  >  S  X.  for  any  h  =  0,1,...,M-1. 
j=l   J    i=l   L 

Further,  let  the  path  in  Bk  first  reach  the  line  Z2  =  N  +  h  (see 

Figure  1)  for  any  h  =  0,1,..., M-l  at  some  node  (w,Nfh)  where  M+k  <  w  <  L. 

Since  the  path  reaches  Z9  =  N  +  h  through  a  vertical  step,  the  (Nfh+l)th 

step  is  just  beginning  at  node  (w,N+h)  while  the  (w+l)th  horizontal  step 

has  already  begun.   Since  the  upper  boundary  must  be  reached  before  the 

right  boundary,  it  is  necessary  that 

L  L 

£X.>  E       Y.    for  M  +  k  <  w  <   L. 

i=w+l    X        j=N+h+l   J 
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This  in  turn  implies  that 

L  L 

(2)     S   X.  >    T.       Y  for  any  h  =  0,1,...,M-1. 

i=M+k+l  1        j=N+h+l  J 

Now  if  the  inequalities  (1)  and  (2)  above  are  added  together,  it  can  be 

k     h 
seen  that  T   >  T   f or  h  =  0, 1, . . . ,M-1 . 
is      D 

c.  If  the  path  taken  was  in  C ,  T  is  larger  than  the  other 
feasible  values.   This  is  proven  analogously  to  part  a. 

d.  If  the  path  taken  was  in  D,  ,  for  some  h  =  0,1,..., M-l,  T 

is  larger  than  the  other  feasible  values.   This  is  proven  analogously 

to  part  b. 

This  concludes  the  proof  of  Theorem  2. 

By  way  of  illustrating  the  result  of  Theorem  2,  an  example  where 

constant  individual  replenishment  times  are  involved  is  discussed. 

It  is  assumed  that  X.  =  A  and  Y.  =  B  for  i,i  =  1,2,...,L  for  the 

i  J 

general  case  of  L  combatant  ships,  distributed  into  queues  of  length 

M  and  N,  associated  with  the  first  and  second  replenishment  ships, 

respectively.   If  neither  N  =  0  nor  M  =  0,  the  following  are  the 

feasible  values  for  T: 

T  =  LA 
A 

T_k  =  (N-k)A  +  (k+l)B  for  k  =  0,1,. . . ,N-1 
Tc  =  LB 

T  h  =  (h+l)A  +  (M-h)B  for  h  =  0,1,..., M-l. 
T,  the  maximum  of  all  feasible  values,  is  clearly  equal  to  either  LA 
or  LB,  depending  upon  the  values  of  A  and  B.   If  A  >  B,  then  T  =  LA;  . 
if  B  >  A,  then  T  =  LB .   If  N  =  0,  all  feasible  paths  are  in  D  and  the 
following  are  the  feasible  values  of  T: 

T  h  =  (h+l)A  +  (L-h)B  for  h  =  0,1,..., L-l. 
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T  is  clearly  equal  to  either  (A  +  LB)  or  (IA  +  B) ,  depending  upon  the 
values  of  A  and  B.   If  A  >  B,  then  T  =  (IA  +  B) ;  if  B  >  A,  then 
T  =  (A  +  LB).   If  M  =  0,  the  preceding  argument  holds  if  the  set  D 
is  replaced  by  the  set  B. 
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IV.   IDENTIFICATION  OF  THE  GENERAL  PROBLEM 

In  the  preceding  section  it  has  been  shown  that  the  actual  total 
replenishment  time,  denoted  T,  is  the  maximum  of  all  feasible  values. 
Thus,  while  the  exact  composition  of  T  as  expressed  in  terms  of  the 
individual  replenishment  times  is  known,  the  general  form  of  its  dis- 
tribution function  as  expressed  in  terms  of  those  for  the  individual 
replenishment  times  remains  unidentified.   In  regard  to  this  last 
statement,  in  this  or  any  succeeding  section  the  following  notation 
has  been  used:   the  distribution  function  of  X.  is  denoted  F. ,  while 

1  1 

that  of  the  Y.  is  denoted  G.  :  if  densities  are  assumed,  they  will  be 
J  J 

denoted  by  the  appropriate  small  letter  with  subscript. 

The  composition  of  T  and  the  distribution  function  of  T  when  L  <  2, 
while  of  little  practical  interest,  are  presented  below  as  an  introduction 
to  the  analytical  techniques  used  in  the  attempt  to  identify  the  general 
distribution  function  of  T,  when  L  >  3: 

a.  if  L  =  1  with  M  =  1  and  N  =  0,  then  T  =  X  +  Y  and  its  dis- 
tribution function  is  the  convolution  of  F,  and  G,  ; 

b.  if  L  =  2  with  M  =  2  and  N  =  0,  then  T  =  X±   +  max  (X^Y^  +  Y2 
and  its  distribution  function  is  the  convolution  of  F..  ,  H_ .  ,  and  G„ , 
where  H9,  denotes  the  distribution  function  of  max  (X_,Y,),  i.e., 

F2G1  i 

c.  if  L  =2  with  M  =  1  and  N  =  1,  then  T  =  max  (X  ,Y.)  +  max  (X2,Y2) 

and  its  distribution  function  is  the  convolution  of  H, .  and  H_„ ,  where 

H. .  denotes  the  distribution  function  of  max  (X.,Y.),  i.e.,  F.G.  for 
n  l   l  '         l  l 

i  =  1,2. 
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It  should  be  noted  that  for  all  cases  investigated,  either  above 

or  beyond  this  point,  an  assumption  that  M  >  N  has  been  made.   The 

symmetrical  nature  of  the  model  has  insured  that  any  results  obtained 

here  are  valid  upon  the  interchange  of  M  with  N,  X.  with  Y.,  and  F. 

1       j       1 

with  G.  for  the  case  when  M  <  N. 
J 

The  general  distribution  function  of  T  can  be  represented  by  the 
following  probability  statements: 

a.  if  N  =  0,  then 

Jl+1        L  . 

P(T  <  t)  =  P  (  Z  X.  +   S  Y.  <  t  for  h  =  0,1,. ..,L-l); 
i=l  L   j=h+l  J  y 

b.  if  N  >  1,  then 

r   L  L        k+1 

P(T  <  t )  =  P  (  £  X .  <  t ;    S   X.+  £  <  t  f or  k  =  0 , 1 , . . . ,N- 1 ; 
i=l  X  "     i=Mfk+l  L   j  =  l  " 

L        h+1         L  . 

EY.<t;  £  X.  +    £   <  t  for  h  =  0, 1, . . . ,M-l). 
j=l  J  '     i=l  L   j=N+h+l  y 

Viewing  the  summation  terms  involved  in  the  two  probability  state- 
ments as  random  variables,  denoted  as  T  ,  it  can  be  noted  that  when 

n 

N  =  0  L  such  random  variables  are  present,  while  if  N  >  1  L  +  2  such 
random  variables  are  present.   Since  it  was  desired  to  use  the  technique 
of  transformation  of  variables  to  simplify  these  probability  statements, 
it  was  obvious  that  additional  random  variables  had  to  be  introduced  in 
order  that  the  one  for  one  transformation  demanded  by  the  technique 

could  be  attained.   Since  the  transformation  was  to  be  from  the  X.  and 

l 

the  Y.  into  the  T  ,  L  artificial  random  variables  would  have  to  be 
J  n' 

introduced  if  N  =  0,  while  L  -  2  would  suffice  if  N  >  1.   These 
artificial  random  variables  were  selected  to  fulfill  two  purposes: 
primarily,  they  insure  a  one-to-one  transformation,  while  secondarily, 
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they  have  no  effect  on  the  appropriate  probability  statements.   For 

L 
the  case  of  N  =  0,  random  variables  of  the  form  £  X.  for  m  =  1,2,...,L 

i=m  L 

are  introduced.   These  summations  are  specifically  less  than  E  X.  +  Y, . 

L        i=l 
For  the  case  of  N  >  1  random  variables  of  the  form  £  X.  for 

L  i=m 

m  =  2,3,...,  M  and  £  Y .  for  n  =  2,3,...,  N  are  introduced .   These 

j=n   J  L  L 


summations   are   specifically    less    than     £  X.    and     £  Y.    respectively. 

i=l   L  j=l  J 

Thus    for   either   case    the    introduction  of   these    particular   summations 

into   the   appropriate    probability  statement   have   no  effect. 

The  whole   set    of  new  random  variables    is   defined  below: 

a.      if  N  =   0, 


T     -     £  X.    for  m  =    1,2, 
m         .         l 
i=m 


h+1  L 

Ti*h+i=  .Vk+  ■  Lyj  forhsa0>1>--->L-1; 

1=1  j=h+l   J 


if  N  >   1, 


T     =     £  X.    for  m  =    1,2,. . . ,M 
m        .         l 

i=m 


L  k+1 

T^i^.i    =  2       X.   +     £  Y.    for  k  =  0,1,... N-l 

>H-k+l         .        .     .    l         .    ,    J 
i=Mfk+l  J  =  l 


TT,      =     £  Y.    for  n   =   1,2,.. . ,N 

h+1  L 

TL+W.h+1  =     S  X.   +  £       Y      for   h   =   0,1,...,M-1. 

1=1  j=N+h+l 

The   Jacobian  matrices    of   these    transformations   are    shown    in  Figures   2, 

3   and   4.      The  value    of   their  determinants    is    the    subject   of   Lemma    1. 

Lemma    1 

The  Jacobian  determinant  of  the  transformation  matrices  shown  in  Figures 

2,  3  and  4  is  equal  to  1,  unless  L  =  2  with  M  =  1  and  N  =  1,  i.e.,  if 
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JACOBIAN  MATRIX   OF   T   TRANSFORMATIONS  WHEN  N  =   0 
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FIGURE   2 
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JACOBIAN  MATRIX  OF  T  TRANSFORMATIONS  WHEN  N  =  1 
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|j|  denotes  the  determinant  of  the  Jacobian  transformation  matrix  J, 
then  |j|  =  1. 
Proof 

a.  Rows  1  through  L  +  N  have  all  O's  to  the  left  of  the  main 
diagonal  elements  and  the  main  diagonal  elements  are  equal  to  1. 

b.  Rows  L  +  N  +  1  through  2L-1  can  be  reduced  to  similar  form  by 
the  application  of  the  formula 

•'      R'T        m  j     •    =  Rt        m  j     •    -   Ri   +  r-j.i    for   i  =   1,2,...,M-1 
L+N+i  L+N+i  1  l+l 

where  R.  denotes  the  ith  row  before  application  of  the  formula,  and 

R'  .  afterwards . 

l 

c.  Row  2L  can  be  reduced  to  a  similar  form  by  application  of  the 
formula 

R'2L=R2L-  Rl  «■-». 

R,2L  =  R2L  -  Rl  +  "Wl  "  RW1  +  RM  "  Rl  +  R2  if  N  "  !•  °r 
R,2L  =  R2L  "  Rl  +  "Wl  "  RL+1  +  RL+2  "  «  >  2. 
The  case  for  L  =  2  with  M  =  1  and  N  =  1  was  specifically  excluded 
from  Lemma  1  since  the  determinant  of  its  Jacobian  transformation 
matrix  equals  zero.   This  results  from  the  inherent  dependence  of  the 
involved  summations.   It  is  noted  that  this  is  the  only  case  in  which 
no  artificial  random  variables  can  be  introduced.   Of  course,  the 
distribution  function  of  T  in  this  case  was  given  earlier  in  this 
section. 

A  second  transformation  of  variables  was  made  at  this  point  to  gain 
further  simplification.   The  following  random  variables  were  defined: 

V  =  T   -  T  ,.   for  m  =  1,2,...,2L-1 
m    m    m+1 

V  =  T 
2L    2L 
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The  determinant  of  the  resulting  Jacobian  transformation  matrix  equals 

1  since  the  matrix  is  triangular  with  all  main  diagonal  elements  equal 

to  1,  i.e.,  all  elements  below  the  main  diagonal  are  O's. 

At  this  point  it  was  decided  to  assume  that  the  distribution 

functions  of  the  individual  replenishment  times  were  continuous  and 

possessed  non-vanishing  density  functions  when  the  X.  or  Y.  are  greater 

than  0.   The  joint  density  function  of  the  X.  and  the  Y.  can  thus  be 
j  j  l  j 

written 

^L        L 

tt  fi(^i)   tt  g.  (y.)  x  >  0  and  y  >  0 
i=l      j=l  J   J  J 

f(x1,.  .  .  ,xL,y][,  .  .  .  ,yL)  =  <)  for  i  ,  j=l  ,2  , .  .  .  ,L 


V. 


otherwise , 


since  it  was  previously  assumed  that  the  X.  and  the  Y.  were  independent 

random  variables. 

Since  the  appropriate  Jacobian  determinants  are  equal  to  1,  the 

ioint  density  function  of  the  V  can  also  be  written  in  the  same 
J  J  n 

manner 


tt  f.(x.)  tt  g.(y.)  x.  >  0  and  y .  >  0 
i=l  L   1  j=l  J   J     1  J 

g(v1,v2,.  .  .  ,v2L)  =|  for  i,j=l,2,...,L 

0  otherwise . 

The  region  in  V-space  for  which  g(v  ,v  , . . . ,v9  )  does  not  vanish  can  be 

represented  in  terms  of  the  V   if  the  transformations  relating  the  X. 
r  m  l 

and  the  Y.  to  the  V  are  used.   These  are  given  below: 
J         m 


if  N  =  0, 


X.  =  V.   for  i  =  1,2,.. . ,L-1 
l    i 

2L 

X.  =  E  V. 

L      -_T   1 

i=L 
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Yj  "Vl  +  V»j     f°r  j   =  W.....L-2 


2L 


b.      if  N  =   1, 


L-l         .    _      1  2L-1 

i=L 


2L-1 
Y     =   -     £       V. 
L  1-1        l 


X.    =  V.      for   i   =   1,2,. .. ,M-1 
l  l 


*M  =  Vl+ 


M  L+l 


2L 

X.    =   -  V.   +  VT   +        £     V. 
L  1  L        .        _    i 

i=L+2 


Y.   -  V,  +  VTJ_.      for   j    =   l,2,...,L-2 
J  J  L+j 

Y  =  V     +  V     +  V  +V 

Vl  1  M  L+l  2L-1 

L-l 

Y.    =   -   2V.    -      S     V.    -   V...   +  VOT 
L  1         .-i  L+l  2L 

i=2 


c.      if  N  >  2 


X.    =  V.       for    i   =   1,2,. . . ,M-1 
l  l 

^"Vk+Wl     fork-  0,1,..., N-2 

XL-1=V1  +  VL-1  +  VL+N 

2L 
XT    -   -  V.   +  VT   +  S        V. 

i=L+N+l 

Yj    =  VL+j      f°r    j    =    1'2"-"N-1 

YN*h   "  Vl  +  VLHHh      for   h  =   0,l,...,M-2 

Y  =  V     +  V  +  V 

L-l  M  L+l  2L-1 

M 

YT    =   -      £     V.    -   V..1   +  V.T 
L  .    .      l  L+l  2L 

i=l 

The   above   statements    thus   allow   the    original   probability   state- 
ments   for   the   distribution   function   of  T   to  be   represented    in   terms 
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of  the  V   in  the  following  form: 
m 

2L 


'(T  <  t)  =  P  T  E  V.  < 


t  for  m  =  1,2 


i=m 


,...,2l) 


provided  that  the  joint  density  function  for  the  V  does  not  vanish 


m 


As   an  example,    the    joint   density   function  of   the  V    ,    expressed  wholly 

in  terms    of   the   V    ,    is    presented   below   for   the   case   N  >  2: 
m  r  — 


M-l 


L-2 


TT    f.(v.)    17    f  .  (V.+VXL^ .      ,)f        ..(v.+V        !+VT,,T) 

.    ,ii.    ..   ix      N+i+1      L-l      l      L-l      L+N 
i=l  i=M 


8(v1}v2 


2L  N-l 


L-2 


'V2L}   = 


V'W.    T^yi).TT18j(VL+j).Vj(Vl-N+VL+j> 

i=L+N+l      j  =  l    J  J    j=N   J      J  J 


M 

8L-l(VVL+l+V2L-l)8L(".E1Vi-VL+l+V2L) 

1=1 


if  condition    (1)   below  exists 


^0  otherwise 

Condition    (1)    is   expressed  below: 


v.    >  0  i   =   1, . . . ,M-1 

l 


vi  +  vtH-i+i>0  l  =  "•••• L"2 


v     +  v  +  v  >  0 

Vl  L-l  L+N 


2L 

-v     +  v     +         S       v.    >  0 
1  L        i=L+N+l    X 


vL+.   >0     j   =   1.....H-1 

Vj+l-N+VL+j   >0      J   =N,...,L-2 
VM+  VL+1  +  V2L-1  >0 


M 

•  S  v 
i=l 


•        "      VT,1       +     VOT       >      0 

i  L+l  2L 


It  is  evident  that  the  above  discussion  has  not  concluded  in  a 
truly  desirable  result,  namely  a  general  form  of  the  distribution 
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function  of  T  that  would  facilitate  ease  of  computation.   The  resulting 

forms,  those  containing  the  T  and  the  V  ,  do,  however,  exhibit  prop- 
'  n  m  f      t- 

erties  that  may  prove  useful  in  future  attempts  in  this  area.   As  the 
matter  now  stands,  not  even  a  relatively  simple  problem,  such  as  that 
when  L  =  3  with  M  =  2  and  N  =  1,  can  be  directly  resolved.   Appendix  A 
depicts  this  particular  problem  as  an  illustration  of  the  general  dif- 
ficulties encountered.   In  retrospect,  it  can  be  stated  that  the 

dependence  relationships  between  the  T  ,  or  those  between  the  V  , 
r  r  n  m 

impose  restrictions  that  cannot  be  overcome  without  some  convenient 
manner  of  handling  the  joint  distribution  function  of  dependent  random 
variables,  such  as  that  developed  in  this  section. 
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V.   A  SIMULATION  OF  THE  UNDERWAY  REPLENISHMENT  OPERATION 

The  preceding  section,  in  which  analytical  methods  were  used, 
failed  to  produce  the  distribution  function  of  T  in  a  form  feasible 
for  specific  computations.   It  was  therefore  decided  to  study  the  dis- 
tribution of  T  by  computer  generation  of  T  realizations  for  a  specific 
case  of  underway  replenishment.   The  primary  intent  of  such  a  study  was 
to  gain  a  better  understanding  of  the  form  of  the  distribution  function 
of  T  for  typical  underway  replenishment  operations. 

The  general  method  used  made  direct  use  of  the  results  of  Theorem 
2,  i.e.,  that  T  is  the  maximum  of  all  feasible  values.   The  underway 
replenishment  operation  was  not  directly  simulated.   Instead,  the 
individual  replenishment  times  were  simulated  and  the  resultant  values 
of  T  calculated,  using  the  results  of  Theorem  2. 

For  the  specific  underway  replenishment  operation,  it  was  decided 
to  use  an  attack  carrier  task  force  replenishment  operation.   The  task 
force  was  selected  to  be  composed  of  eight  combatant  ships,  an  attack 
aircraft  carrier  (CVA) ,  a  light  guided  missile  cruiser  (CLG) ,  two 
guided  missile  frigates  (DLG) ,  and  four  destroyers  (DD) .   The  two  ship 
replenishment  group  was  selected  to  be  composed  of  a  fleet  oiler  (AO) 
and  an  ammunition  ship  (AE) .   This  particular  selection  was  made  because 
approximations  for  the  distribution  functions  of  the  individual  replenish- 
ment times  for  the  above  type  ships  were  presented  by  Besecker  in  [3]. 
Thus  Erlang  distribution  functions  were  assumed  for  the  individual 
replenishment  times;  the  specific  gamma  parameters,  denoted  k  and  \, 
are  shown  in  Table  I.   The  general  form  of  the  Erlang  distribution 
function  is  given  in  Appendix  B. 
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PARAMETERS  FOR  ERLANG  DISTRIBUTION  FUNCTIONS 

FOR 
INDIVIDUAL  REPLENISHMENT  TIMES 

REPLENISHMENT  SHIPS 

AO 

AE 

k 

X 

k 

X 

C 
0 

M  S 
B  H 
A   I 
T   P 
A  S 
N 
T 

CVA 
CLG 
DLG 
DD 

11 
7 
7 
4 

0.09 
0.11 
0.11 
0.05 

3 
2 
2 
3 

0.03 
0.02 
0.02 
0.06 

TABLE  I 

In  accordance  with  the  model,  the  AO  and  AE  were  designated  the 

first  and  second  replenishment  ships,  respectively.   The  AO  is  then 

associated  with  the  X.,  while  the  AE  is  associated  with  the  Y..   It 

J-  J 

was  decided  to  establish  a  constant  replenishment  sequence  at  the  AO. 
This  sequence  is  CVA,  DLG,  DD,  DD,  CLG,  DLG,  DD,  and  DD.   The  replenish- 
ment sequence  at  the  AE  would  depend  upon  the  M  and  N  combination  being 
used  for  the  generation  of  a  T  realization.   The  feasible  sequences 
are  shown  in  Table  II. 

The  individual  replenishment  times  were  simulated  by  the  inverse 
transformation  method  discussed  in  Appendix  B.  The  property  that  an 
Erlang  random  variable  is  the  sum  of  k  random  variables,  each  expon- 
entially distributed  with  parameter  X,    facilitates  this  method.   To 
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insure  that  the  random  variates  generated  were  from  the  assigned  dis- 
tributions, chi-square  goodness  of  fit  tests  were  made  at  the  0.01 
level  of  significance.   This  method  is  discussed  in  Appendix  B.   For 
all  16  distributions  of  the  individual  replenishment  times,  the  null 
hypotheses  that  the  generated  random  variates  came  from  the  appropriate 
distributions  were  of  course  accepted. 


- 

REPLENISHMENT 

SEQUENCE  AT 

AE 

M 

8 

7 

6 

5 

4 

3 

2 

1 

0 

N 

0 

1 

2 

3 

4 

5 

6 

7 

8 

R 

CVA 

DD 

DD 

DLG 

CLG 

DD 

DD 

DLG 

CVA 

E 
P     S 

DLG 

CVA 

DD 

DD 

DLG 

CLG 

DD 

DD 

DLG 

L     E 
E      Q 

DD 

DLG 

CVA 

DD 

DD 

DLG 

CLG 

DD 

DD 

N     U 
I      E 

DD 

DD 

DLG 

CVA 

DD 

DD 

DLG 

CLG 

DD 

S      N 
H     C 

CLG 

DD 

DD 

DLG 

CVA 

DD 

DD 

DLG 

CLG 

M      E 
E 

DLG 

CLG 

DD 

DD 

DLG 

CVA 

DD 

DD 

DLG 

N 
T 

DD 

DLG 

CLG 

DD 

DD 

DLG 

CVA 

DD 

DD 

DD 

DD 

DLG 

CLG 

DD 

DD 

DLG 

CVA 

DD 

TABLE  II 

A  special  feature  was  provided  in  the  implementation  of  the  com- 
puter program  to  enable  study  of  the  dependence  of  T  on  the  M  and  N 
combination  being  used.  The  program  provided  that  upon  the  nth  gen- 
eration of  a  T  realization  for  a  particular  M  and  N  combination,  all 
ships  would  have  the  same  individual  replenishment  times  at  both  re- 
plenishment ships  as  for  any  other  M  and  N  combination,  i.e.,  on  the 
50th  generation  of  a  T  realization  for  M  =  3  and  N  =  5,  the  CLG  would 
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have  the  same  replenishment  times  at  both  the  AO  and  AE  as  for  the  50th 
generation  of  a  T  realization  for  any  other  M  and  N  combination  of  the 
combatant  ships. 

Five  hundred  T  realizations  were  generated  for  each  of  the  nine 
feasible  M  and  N  combinations.   The  random  walk  for  a  typical  T 
realization  is  shown  in  Figure  5  as  a  random  walk  within  the  constraint 
boundary  of  the  L  x  L  grid  and  in  Figure  6  as  a  random  walk  in  two- 
dimensional  time  space.   The  T  realization  obtained  by  performing  the 
random  walk  agrees  with  that  obtained  through  an  application  of 
Theorem  2 . 

For  each  M  and  N  combination,  the  generated  data  were  used  to  con- 
struct frequency  histograms  of  T  realizations.   Figure  7  is  a  typical 
frequency  histogram.   From  these  histograms,  one  basic  observation 
could  be  made:   the  frequency  histograms  showed  definite  positive 
skewness ,  i.e.,  they  were  skewed  to  the  right.   The  obvious  contention 
at  this  point  was  that  the  distribution  function  had  the  same  character- 
istic shape  as  did  those  of  the  individual  replenishment  times,  i.e., 
the  shape  of  a  gamma  distribution  function. 

A  second  observation  could  be  made  once  the  T  sample  means  for  the 
various  M  and  N  combinations  were  calculated.   The  results  are  similar 
to  those  cited  by  Milch  and  Waggoner  in  [4]  for  the  expected  values  of 
T  when  the  X.  and  Y.  are  exponentially  distributed  with  parameters  X 
and  [i,    respectively.   The  minimal  T  sample  mean  occurred  when  M  and  N 
were  equal  and  a  range  of  M  and  N  combinations  exists  in  which  there 
is  only  a  slight  variation  of  the  T  sample  mean  values.   The  T  sample 
means  are  plotted  against  M  and  N  combinations  in  Figure  8. 

It  was  decided  to  investigate  the  contention  that  the  distribution 
function  of  T  was  from  the  family  of  gamma  distributions.   Estimates  of 
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FIGURE  5 
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A       A 

necessary  gamma  distribution  parameters  k  and  \,  denoted  k  and  \,   were 
obtained  by  the  method  of  moments  as  shown  in  Appendix  B.   The  basic 
null  hypothesis  to  be  tested  would  be  that  the  T  realizations  for  a 
particular  M  and  N  combination  came  from  a  gamma  distribution  with 

A        A 

parameters  equal  to  k  and  \.   The  chi-square  method  was  choosen  for 
testing  goodness  of  fit.   Despite  the  generally  acceptable  results, 
eight  out  of  nine  null  hypotheses  being  acceptable  at  the  0.01  level 
of  significance,  it  was  evident  that  better  results  could  be  achieved. 
By  studying  superimposed  frequency  histograms  of  the  observed  T 
realizations  and  the  expected  T  realizations  under  the  gamma  distri- 
bution of  the  null  hypotheses  for  both  the  accepted  and  rejected  cases, 
two  observations  could  be  made: 

a.  the  peak  ordinate  value  for  the  number  of  observed  T  realizations 
was  higher  than  that  of  the  peak  ordinate  value  for  the  number  of 
expected  T  realizations  for  all  M  and  N  combinations  ;  and 

b.  the  mode,  the  absissca  for  the  peak  ordinate  value,  for  the 
number  of  observed  T  realizations  either  occurred  to  the  right  of  or 
approximately  coincided  with  the  mode  for  the  number  of  expected  T 
realizations  for  all  M  and  N  combinations. 

A  typical  superimposed  frequency  histogram  is  shown  in  Figure  9. 

From  these  observations,  the  following  generalization  could  be 
made  : 

the  true  density  function  of  T  has  a  maximum  value  that  occurs 
above  and  to  the  right  of  a  gamma  density  function  with  parameters 

A        A 

k  and  X . 

Thus  it  was  decided  to  use  a  method  of  estimating  gamma  parameters 
that  utilized  the  maximum  value  of  the  gamma  density  function.   This 
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method,  termed  the  peak  value  method,  is  developed  in  Appendix  B. 
Estimates  of  k  and  X,  as  determined  by  this  method,  are  denoted  k*  and 
X*,  respectively.   Success  was  readily  achieved  by  this  method.   The 
basic  null  hypothesis,  that  T  for  a  particular  M  and  N  combination  was 
distributed  by  a  gamma  distribution  with  parameters  k*  and  X*,  was 
acceptable  at  the  0.01  level  of  significance  by  the  chi-square  method 
for  all  nine  M  and  N  combinations.   The  results,  including  those 
obtained  by  the  method  of  moments,  are  shown  in  Table  III.   It  could 
be  concluded  that  the  distribution  function  of  T  could  be  successfully 
approximated  by  a  gamma  distribution  function  with  parameters  k*  and  X*. 


RESULTS    OF  METHOD   OF  MOMENTS   AND   PEAK  VALUE 

ESTIMATION 

M/N 

T 

sT2 

A 

k 

A 

X 

k* 

X* 

8/0  # 

13.12 

4.24 

40.58 

3.09 

46.19 

3.52 

7/1 

12.28 

4.08 

36.95 

3.01 

41.32 

3.36 

6/2 

11.83 

3.76 

37.24 

3.15 

44.63 

3.77 

5/3 

11.52 

3.18 

41.76 

3.63 

48.66 

4.23 

4/4 

11.49 

3.17 

41.67 

3.63 

47.70 

4.15 

3/5 

11.50 

3.19 

41.46 

3.60 

46.34 

4.03 

2/6 

11.53 

3.22 

41.28 

3.58 

46.60 

4.04 

1/7 

11.95 

3.97 

35.91 

3.00 

45.06 

3.77 

0/8 

13.46 

4.71 

38.43 

2.86 

46.93 

3.49 

#:   M  and  N  combination  that  resulted  in  rejection  of  null 
hypothesis  under  method  of  moments  testing 

TABLE  III 
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APPENDIX  A 

ILLUSTRATION   OF    PROBLEM  WHEN   L  =  3   WITH  M  =  2   AND   N  =    1 
The  elements    of   this    particular   problem  are    presented  below   in   the 
context    of   Section    IV: 

a.      Considering   the    feasible    paths,    this    problem   involves    random 
variables    of   the    form: 

TA=X1  +  X2  +  X3 

X3  +  Yl 

Yl  +  Y2  +  Y3 


«S 

= 

Tc 

= 

i 

= 

xi 

i 

= 

xi 

+ 

X2 

+  Y2  +  V3 

+  Y3- 

The  distribution   function   of   T  can  be  written 

P(T  <    t)    =   P(X  +X2+X3   <   t,    X3+YL  <   t,   Y1+Y2+Y3   <   t, 

Xl+Y2+Y3  -  fc'   Xl+X2+Y3  -  t)' 
b.      After    the   addition   of   the   appropriate   artificial   random  variable, 
the    problem   involves    random  variables    of   the   form: 

i.  1      =    X,     T    X.,j     T    X« 

T2   =  X2  +  X3 

T3=  X3+Yl 

T,   =  Yx  +  Y2  +  Y3 

T5   =  Xl  +  Y2  +  Y3 

T6=X1+X2  +Y3. 
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The  distribution  function  of  T  can  be  written 

P(T  <  t)  =  P(T.  <  t  for  i  =  1,2,. ..,6). 

c.   After  the  transformation  from  the  T  to  the  V  ,  the  problem 

n         m 

involves  random  variables  of  the  form: 


X, 


V3  = 


-Y   -Y 
2   3 


VA  =-Xl 


+Y, 


v5=  -x2 


+Y, 


v6  =  x1+x2 


+Y3- 


The  distribution  function  of  T  can  be  written 


P(T  <  t)  =  P  T.   V   t  for  i  =  1,2,. ..,6). 
—  .  m 

m=i 

d.   Considering  the  inverse  transformations  from  the  V  to  the  X. 

mi 


and    the  Y.,    the    problem   involves    the    forms 

xi=  v 


x2=  V 


x3  =  -V 


v1=    V 


Y  =  V 


Y3  =-2V 


+  V, 


+  v. 


+  V, 


+  Vc  +  V, 

J      D 


+  v, 


+  vr 


-  v, 


+  v4  +  v5 


-  v, 


+  v 


e.   The  joint  density  function  of  the  V  ,  as  wholly  expressed  in 

terms  of  the  V  ,  can  be  written 
m 


A7 


•f  x  (v  l )  f  2  (v  L+v2+v4  )  f  3  (  - v  1+v3+v5+v6  )   x 

g(v1,v2,v3,v4,v5,v6)    =   JMVVMVVVV^^VVVV 

if   condition    (1)    below  holds 
0  otherwise 

Condition    (1)    is   as    follows 

V,  >  0 


VV2  +  V4  >  ° 


V     +  V     +  V     >  V 

v3   -t-  v5  f   v6   >  V]L 


vx  +  v4  >  0 


Vl  +  V2  +  V4  +  V5  >  ° 


v6  >  2Vi  +  v2         +  v4 
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APPENDIX  B 

RELEVANT  THEORETICAL  AND  STATISTICAL  FORMS 
1 .   The  Gamma  Distribution 

The  gamma  distribution  function  is  given  by 


F(x)  = 


f  ,k  k-1   -\x  . 

*   x    e     dx    c  . 
sry.   . for  x  >  0  and  k  >  0, 

r(k) 


/k-1  -x  2 

x    e  ""  dx.   Let  (j,,  m,  and  a  denote  the  mean,  the  mode, 

o 
and  the  variance  of  the  gamma  distribution,  then 

M  =  k/X ; 

m  =  (k-1) A ;  and 

a  =  kA  • 

If  k  is  a  positive  integer,  then  the  resulting  gamma  distribution  is 

also  termed  an  Erlang  distribution.   The  Erlang  distribution  can  be 

shown  to  be  that  for  the  sum  of  k  identically  distributed  independent 

random  variables  from  an  exponential  distribution  with  parameter  X. 

If  k  is  of  the  form  n/2  with  n  a  positive  integer  and  X   equals  1/2, 

then  the  resulting  gamma  distribution  is  also  termed  a  chi-square 

distribution. 

2 .   Determination  of  Erlang  Random  Variables  by  the  Inverse  Trans- 
formation Method 

Let  F(x)  denote  the  distribution  function  for  a  random  variable  X 

and  let  r  denote  a  random  number  from  the  uniform  distribution  on  the 

interval  0  to  1.   If  r  =  F(x),  then  x  =  F(r)   ,  where  F(r)   denotes 

the  inverse  transformation  of  F(x).   If  X.  is  exponentially  distributed 

The  transformation  of  a  random  variable  X  into  a  uniformly  distributed 
random  variable  r  =  F(x)  is  a  well  known  result  in  probability  theory. 
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with  parameter  X,  then 

F(x.)  =  l-e"^Xi;  or  if  F(x.)  =  r., 
iy  11 

then  r.  =  1-e   1; 

1 

-Xx 
thus  e    i  =  1  -  r.  :  but  if  r.  is  a  realization  from  the  uniform  dis- 

l  l 

tribution  on  the  interval  0  to  1,  then  1  -  r.  =  r. *  is  also  a  realization 

l    i 

on  the  same  interval.   Thus  this  may  be  rewritten 

-Xx 
e    i  =  r. '  :  solving  now  for  x. , 
l  i 

x.  =  -  1/X  Inr.'. 
l  l 

k 

If  X  =  E  X. ,  then  X  has  an  Erlang  distribution  with  parameters  k  and 
i=l   X 


X ;  thus 


X  =  -  1/X   E   In  r . ' ;  which  may  be  rewritten 
i=l     X 


X  =  -  1/X  In  C  tt  r.  *\  where 


k 

TT 

1=1 


F(r)~   =  -  1/X  In  (    TT  r.  M  with  r  =  tt  r.  '  . 

i=l  i=l 

3.   Test  for  Goodness  of  Fit  by  Chi-square  Criteria 

Let  x.  be  a  realization  from  the  sampled  population;  let  0.  denote 

the  number  of  realizations  that  fall  into  class  j  for  j  =  l,2,...,n, 

where  n  is  equal  to  some  arbitrary  number  of  subdivisions  of  the  range 

of  the  assumed  distribution  of  X.;  let  E.  denote  the  expected  number 

i       J 

of  realizations  in  class  j  under  the  assumed  distribution;  then  the 
statistic 

n   (E.  -  0.)2 

I\s1-LirL- 

has  an  approximate  chi-squared  distribution  with  parameter  (n-m-l)/2, 
where  m  denotes  the  number  of  assumed  distribution  parameters  estimated 
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with  sample  statistics,   n-m-1  is  called  the  degree  of  freedom  associ- 
ated with  X;  for  a  choosen  significance  level  the  value  of  X  can  be 
compared  with  the  appropriate  table  value  such  as  found  in  Appendix  4 
of  [5];  if  X  is  greater  than  the  table  value,  then  the  null  hypotheses 
that  the  x.  came  from  the  assumed  distribution  is  rejected;  it  is 
accepted  otherwise. 

4.  Estimation  of  Gamma  Parameters  by  the  Method  of  Moments 

Let  T  be  a  random  variable  from  a  gamma  distribution  with  parameters 

k  and  X.      Let  |_l  denote  the  expected  value  of  T,  T  denote  the  T  sample 

2  2 

mean  value,  a  denote  the  variance  of  T,  and  S   denote  the  T  sample 

n     -  2 
variance  as  calculated  by  l/(n-l)  Z  (T.-T)  ,  where  n  is  the  number  of 

i=l   L 

T  sample  values.   Since 

HT  =  k/X ,  and 

2   .,,2 
Ot  =  k/K    , 

A        A 

it  can  be  shown  that  k  and  X,    estimates  of  k  and  A. ,  respectively,  can 

be  calculated  by 

a   —2   o 

k  =  r/SJ  and 

X  =  f/s*. 

2 
if  the  mean  |j,  and  variance  o_  are  respectively  estimated  by  the  sample 

mean  and  variance. 

5 .  Estimation  of  Gamma  Parameters  by  the  Peak  Value  Method 

Let  T  be  a  random  variable  from  a  gamma  distribution  with  parameters 
k  and  X.      Let  u  denote  the  expected  value  of  T,  T  denote  the  T  sample 
mean  value,  p  denote  the  maximum  value,  i.e.,  the  peak,  of  the  relevant 


gamma  density  function,  and  m  denote  the  mode  of  the  distribution  of  T, 


A    A 

i.e.,  the  T  value  at  which  p  occurs.   Let  \±    ,  p,  k*  and  X*   denote 


51 


estimates  of  [i_,  p,  k  and  X,    respectively.   Since  p  occurs  at  m  and 
since  m  =  (k-l)A,  p  can  be  evaluated  as 

(k-pfc-1  e'-k 

p  =         r(k)       • 

Now  by  applying  an  approximation  formula  for  T(k)  as  found  in  [6], 

i.e., 

. —     _k     k-1/2 
F(k)   =  Jlr\  e     '  k  for  k   large   and   not  an   integer, 

p  can  be   approximated   as 


(l-l/k)k-1 


2TTk 

k-1  -1 

Since  for  large  k,  (1-1/k)    can  be  approximated  by  e   ,  p  can  be 

approximated  as 

V2nk 

A  A 

To  estimate  k  and  X,   \i     and  p  will  be  used.   Since  |a  =  k/X   and  since 
p  can  be  approximated  by  formula  (1),  by  combining  these  equations  and 
solving  for  k,  it  is  easily  seen  that 

2  2 

k  ~  2np  m^. 

Now,  knowing  an  approximation  for  k,  an  approximation  of  X    can  be 

obtained  from  X   =  k/^a  .   Now,  by  letting  (^  =  T  and  p  =  max  O./nw, 

where  max  0.  is  the  number  of  T  realizations  in  the  modal  observation 
J 

class  of  the  T  sample  frequency  histogram,  n  is  the  number  of  T 
realizations  in  the  sample,  and  w  is  the  width  of  the  modal  class  in 
the  frequency  histogram,  k*  and  X*   can  then  be  evaluated  as 

k*  =  2irp2T2 

X*   =  k*/T. 
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